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ABSTRACT 

A  study  of  fuel  disintegration  and  mixing  in  a  supercritical  environment  (relative  to  the 
fuel)  was  undertaken  in  order  to  determine  parameter  regimes  advantageous  to  mixing.  The 
adopted  approach  was  based  on  the  development  of  a  model  describing  a  supercritical,  temporal 
mixing  layer  and  on  the  utilization  of  Direct  Numerical  Simulations  (DNS)  as  an  investigation 
method. 

During  the  current  year,  the  database  from  a  DNS  of  a  temporal,  binary  species,  initially 
density  stratified,  excited  mixing  layer  that  did  not  reach  transition  was  enlarged  and  scrutinized 
to  determine  the  reasons  for  the  lack  of  transition.  The  scrutiny  was  based  on  the  development  of 
the  entropy  equation  for  a  supercritical  fluid,  and  on  the  examination  of  the  various  contributions 
to  the  irreversible  entropy  production  (the  dissipation).  The  results  showed  that  the  lack  of 
transition  was  due  to  the  formation  of  regions  of  large  density  gradients  which  damped  the 
emerging  turbulent  scales.  To  obtain  transition,  the  Reynolds  number,  as  well  as  the  amplitude  of 
the  initial  excitation  were  increased,  producing  two  simulations  in  which  transition  was  attained; 
two  new  databases  were  thus  created,  for  initial  Reynolds  numbers  of  500  and  600,  respectively. 
Since  vorticity  is  known  to  play  a  crucial  role  in  turbulence,  the  two  new  databases  were 
examined  to  determine  the  phenomena  responsible  for  vorticity  production.  Examination  of  the 
vorticity  budget  at  transition  revealed  that  most  of  the  spanwise  vorticity  production  was  due  to 
stretching/tilting,  and  that  most  of  the  contribution  to  the  vorticity  magnitude  was  due  to  the 
viscous  terms.  It  was  also  shown  that  the  regions  of  large  density  gradients  created  during  roll-up 
and  pairing  persisted  at  transition,  and  that  the  fluid  in  these  regions  was  non-ideal  and  consisted 
primarily  of  the  entrained  fluid,  with  small  amounts  of  entraining  fluid  dissolved  in  it. 

TECHNICAL  DISCUSSION 

The  first  year  of  this  investigation  was  devoted  to  developing  a  database  of  a  three- 
dimensional,  supercritical,  transitional  temporal  mixing  layer.  The  mixing  layer  was  composed 
of  a  lower  stream  containing  the  fuel  (n-heptane)  and  the  upper  stream  containing  nitrogen; 
supercriticality  was  initially  measured  with  respect  to  the  pure  fuel.  The  mixing  layer  was  chosen 
as  the  simplest  configuration  from  which  one  may  obtain  results  that  will  establish,  through 
DNS,  a  framework  for  understanding  mixing  of  supercritical  fluids  at  small  scales. 

Following  previous  modeling  performed  at  JPL  for  isolated  fluid  drops  (Harstad  and 
Bellan,  1998,  2000),  the  model  of  supercritical  behavior  was  based  upon  fluctuation-dissipation 


theory  (Keizer,  1987).  The  advantage  of  this  theory'  was  that  it  inherently  accounted  for 
nonequilibrium  processes  and  naturally  led  to  the  most  general  fluid  equations  by  relating  the 
partial  molar  fluxes  and  the  heat  flux  to  thermodynamic  quantities.  Thus,  Soret  and  Dufour 
effects  (potentially  important  at  supercritical  conditions)  were  taken  into  account  from  first 
principles  through  the  transport  matrix  where  they  complemented  the  traditional  Pick's  mass 
diffusion  and  Fourier  thermal  diffusion  terms.  Results  from  the  recent  supercritical  simulations 
of  Harstad  and  Bellan,  2000,  using  heptane  drops  in  nitrogen  have  already  been  validated  with 
the  microgravity  data  of  Nomura,  1996,  within  an  average  of  15%,  thereby  showing  that  the  fluid 
model  was  correct.  This  model  was  here  adapted  to  the  shear  layer  configuration.  However,  there 
were  specific  issues  that  were  considered  in  the  DNS  of  a  three-dimensional  (3D)  shear  layer  and 
that  had  not  been  relevant  to  the  drop  model: 

(a)  In  contrast  to  the  drop  study,  the  treatment  of  the  transport  coefficients  cannot  be 
exact  in  the  context  of  DNS  since  the  Reynolds  number.  Re,  would  be  far  too  large  to  make  the 
computations  feasible.  Since  it  was  required  that  the  Batchelor  scales  (the  smallest  scales 
associated  with  scalar  mixing,  tjg  ~  tj^J Sc  where  rjf^  was  the  Kolmogorov  scale  and  Sc  was  the 
Schmidt  number)  be  resolved,  this  dictated  the  maximum  Sc  which  could  be  used  in  the 
particular  calculation.  For  example,  for  heptane/nitrogen.  Sc  was  too  large  for  any  domain  size 
since  heptane  had  a  liquid-like  density.  The  strategy  was  then  to  calculate  these  numbers  based 
upon  the  observations  that  the  viscosity  was  primarily  a  function  of  T,  whereas  Sc  and  Pr  were 
primarily  functions  of  the  mass  fractions,  and  to  find  approximate  fits  for  these  quantities  in  the 
range  of  parametric  interest. 

(b)  It  was  here  unpractical  to  use  the  very  accurate  and  computationally  efficient  non 
ideal  EOS’s  employed  in  Harstad  and  Bellan,  2000,  because  these  could  still  be  too 
computationally  intensive  when  performing  3D  simulations.  Instead,  the  Peng-Robinson  EOS 
was  used,  based  on  pure  species  reference  states  accurate  to  better  than  1%  relative  error  (on  pure 
reference  states)  through  comparisons  with  highly  accurate  EOS’s  (Harstad  et  al.,  1997)  over  the 
range  of  variables  used  in  this  study. 

The  shear  layer  equations  were  coded  using  a  fourth  order  Runge-Kutta  scheme  for 
advancing  in  time,  and  an  eight-order  compact  method  for  spatial  discretization.  The 
implemented  boundary  conditions  were  those  of  Poinsot  and  Lele,  1992,  for  compressible  flows. 
Since  those  boundary  conditions  were  developed  for  perfect  rather  than  real  gases,  they  were 
adapted  for  the  present  real  gas  computations  by  calculating  an  ‘equivalent  gas  constant’  from 
the  EOS;  the  accuracy  of  this  simplistic  substitution  is  discussed  in  more  detail  below.  Also,  the 
layer  was  perturbed  in  the  mean  velocity  using  the  most  unstable  wavelength  of  Moser  and 
Rogers,  1991,  found  in  the  context  of  atmospheric,  incompressible  gases. 

Results  from  the  3D  shear  layer  (initial  conditions:  Reynolds  number  of  400,  Mach 
number  of  0.4,  60  atm,  lower  stream  at  600K,  and  upper  stream  at  lOOOK)  showed  that  the 
density  stratification  was  stabilizing  the  layer  and  preventing  it  from  achieving  the  transitional 
regime  at  a  nondimensional  time  at  which  a  transitional  state  had  been  reached  by  a  layer  whose 
lower  stream  was  initially  laden  with  drops.  One  intriguing  feature  of  the  layer  was  the 
appearance  of  regions  of  high  density  gradient  magnitude  which  were  not  merely  the  distortion 
of  the  original  density  interface.  The  appearance  and  geometry  of  these  regions  of  high  density 
gradient  magnitude  was  intriguing  because  it  corresponded  to  the  wispy  threads  of  fluid 
emanating  from  a  supercritical  nitrogen  jet  in  gaseous  nitrogen  as  observed  by  Chehroudi  et  al., 
1999.  These  features  were  absent  at  subcritical  conditions  either  in  the  experiment  or  in 
atmospheric  air  DNS  simulations,  indicating  that  they  may  be  specific  to  the  supercritical 
situation.  These  results  were  all  documented  in  Miller,  Harstad  and  Bellan,  2000. 


Since  the  layer  did  not  achieve  transition  in  a  physical  time  comparable  to  that,  for 
example,  of  a  drop  laden  mixing  layer,  the  calculations  were  pursued  with  the  goal  of  achieving 
transition  at  a  later  physical  time.  However,  this  did  not  occur  for  the  chosen  initial  conditions, 
and  moreover,  the  layer  reached  a  culminating  point  in  the  averaged  positive  spanwise  vorticity. 
The  averaged  positive  spanwise  vorticity  was  considered  to  be  a  good  indicator  of  transition 
since  the  initial  mean  velocity  profile  was  such  that  the  spanwise  vorticity  was  negative 
everywhere;  transition  would  occur  if  the  positive  spanwise  vorticity  exhibited  and  continued  to 
display  a  large  increasing  rate.  The  results  illustrated  in  Fig.  1,  together  with  the  momentum 
thickness,  the  product  thickness  and  the  enstrophy  (another  indicator  of  transition),  indicated  the 
lack  of  transition. 

To  investigate  the  origin  of  the  stability  preventing  transition  in  a  physical  time  larger 
than  that,  for  example,  of  a  drop  laden  mixing  layer,  the  entropy  equation  was  deri\’ed  for  a 
supercritical  fluid  and  the  irreversible  entropy  (i.e.  the  dissipation)  production  was  examined.  For 
supercritical  fluids,  the  dissipation  was  shown  to  contain  three  contributions:  viscous,  Fourier 

heat  flux  and  molar  flux  contribution. 
The  molar  flux  contribution  contained 
both  the  effect  of  the  mixture 
nonideality  and  that  of  the  Soret  term. 
The  database  from  the  long-time  run 
was  used  to  compute  the  budget  of  the 
small  scale  dissipation  by  calculating 
the  difference  between  the  DNS  (i.e. 
unfiltered)  and  that  of  the  filtered  flow 
field.  The  results  from  plane  averages 
and  RMS  of  the  three  terms  showed 
that  the  dominant  contribution  to  the 
small  scale  average  dissipation  was 
from  the  viscous  term,  but  that  the 
major  RMS  contribution  was  from  the 
molar  flux  term  at  all  times  except  at 
the  culmination  of  the  averaged 
positive  spanwise  vorticity.  Since  the  molar  flux  contribution  to  the  dissipation  contained  six 
terms,  a  budget  of  the  molar  flux  dissipation  was  calculated  to  identify  which  of  the  terms  was 
the  principal  contributor  to  the  molar  flux  dissipation  magnitude.  It  was  found  that  the  non¬ 
ideality  term  was  the  major  contributor  in  the  heptane  stream,  whereas  the  nonideality  and  Soret 
terms  tended  to  be  of  same  order  of  magnitude  in  the  nitrogen  stream.  These  calculations  as  well 
as  contour  plots  showed  that  the  RMS  of  the  molar  flux  dissipation  correlated  with  the  regions  of 
maximum  density  gradient  magnitude.  Moreover,  it  was  also  noticed  that  the  DNS  dissipation 
RMS,  as  well  as  each  of  its  contributions  was  larger  than  the  respective  average  value,  indicating 
that  there  was  considerable  backscatter  in  the  flow. 

The  results  concerning  the  hindering  of  transition  due  to  the  presence  of  a  large  density 
gradient  region  were  found  consistent  with  the  experimental  results  of  Hannoun  et  al.,  1988,  who 
found  that  large  eddies  impinging  on  a  density  interface  ‘bounce’  back  without  significant 
entertainment  of  unstirred  fluid  instead  of  overturning.  Since  the  layer  growth  is  known  to 
depend  primarily  on  entertaimnent,  and  since  this  accelerated  growth  is  known  to  promote  the 
appearance  and  evolution  of  the  small  turbulent  scales,  this  damping  mechanism  may  be 
enhancing  the  effect  of  density  stratification  in  hindering  transition. 


The  calculation  of  the  dissipation  was  valuable  not  only  in  elucidating  the  responsible 
mechanisms  impacting  the  stability  characteristics  of  the  flow,  but  also  in  pointing  out  that  at 
large  simulation  times  the  mass  fractions  and  temperature  displayed  large  gradients  at  the  upper 
and  lower  boundaries  of  the  computational  domain;  these  were  the  only  locations  where  non¬ 
periodic  boundary  conditions  were  applied.  This  provided  the  motivation  to  re-examine  the 
implementation  of  the  boundary  conditions  according  to  Poinsot  and  Lele,  1992.  and  Baum  et  al., 
1994.  Thus,  the  ‘characteristic’  boundary  conditions  were  derived  for  a  general  fluid  having  a 
real  gas  equation  of  state.  This  derivation  was  consistent  with  the  fact  that  for  Rc  oo  the 
conservation  equations  including  the  Soret  and  Dufour  effects  were  an  incomplete  elliptic  set.  In 
fact,  the  elliptic  character  was  determined  to  be  the  result  of  the  value  of  the  effective  thermal 
conductivity  which  was  always  enhanced,  while  the  effective  diffusivity  was  always  reduced  (see 
Harstad  and  Bellan,  1999)  when  comparing  to  atmospheric  flows.  Simulations  of  a  one¬ 
dimensional  propagation  of  acoustic  waves  in  a  two-dimensional  domain  using  subsonic  non¬ 
reflecting  boundaries  were  conducted  with  the  new,  fundamental  boundary  conditions,  and 
compared  with  the  previous  simulations  based  on  the  ‘simplistic’  effective  gas  constant  (see 
above).  The  results  revealed  that  whereas  using  the  fundamental  analysis  induced  the  acoustic 
waves  to  properly  pass  through  the  boundaries  without  reflections,  the  use  of  the  simplistic 
approach  yielded  significant  reflections  at  the  boundaries.  Most  discrepancies  occurred  when  a 
source  term  was  added  to  the  mass  fraction  equations,  where  additional  to  the  reflections  at  the 
boundaries,  the  solution  within  the  entire  domain  was  different  in  the  simplistic  and  fundamental 
approaches.  This  work  will  be  enlarged  in  the  near  future  to  consider  two-dimensional  wave 
propagation,  and  the  results  will  be  further  documented  in  a  manuscript. 

The  fundamental  boundary  conditions  were  implemented  in  the  code,  the  discretization 
was  changed  to  a  sixth-order  compact  scheme,  and  two  simulations  were  performed  for  initial 
Reynolds  numbers  of  500  and  600,  respectively.  Given  the  large  initial  density  stratification 
(12.88  for  the  initial  conditions  of  the  calculation:  60  atm,  lower  stream  at  600K,  upper  stream  at 
lOOOK,  and  Mach  number  of  0.4),  it  was  deemed  important  to  increase  the  amplitude  of  the  layer 
excitation,  so  as  to  increase  the  entertainment  of  the  heavier,  lower  stream  n-heptane.  With  these 
new  initial  and  boundary  conditions,  transition  was  achieved  in  both  simulations.  Portrayed  in 
Fig.  2  are  the  global  characteristics  of  these  two  layers. 


Figure  2  Time  Evolution  of  Global  Quantities  for  Reo=500  and  Reo=600  a)  Momentum  Thickness  b)  Product 
Thickness  c)  Positive  Spanwise  Vorticity  d)  Enstrophy 

The  momentum  thickness  based  Reynolds  numbers  were  1250  and  1452,  respectively,  at 
transition.  These  values  are  smaller  than  the  traditional  magnitudes  of  measured  transitional 
Reynolds  numbers  due  to  the  reference  viscosity  being  here  four  orders  of  magnitude  larger  than 
the  species  viscosity;  this  is  an  artifact  of  DNS  designed  to  allow  the  attainment  of  large 
Reynolds  numbers  in  domains  of  physical  interest. 

Since  vorticity  is  known  to  play  a  major  role  in  turbulent  flows,  the  phenomena 
predominantly  responsible  for  vorticity  production  were  identified  by  analyzing  the  budget  of  the 
vorticity  equation  at  the  transitional  state.  It  was  found  that,  on  average,  the  stretching/tiltin^ 
term  dominated  all  other  effects  in  the  production  of  spanwise  vorticity,  but  that  the  viscous  term 
dominated  the  vorticity  magnitude  production.  In  this  respect,  the  major  difference  identified 
between  the  layer  that  did  not  reach  transition  and  the  present  case  was  the  dominance  of  the 
stretching/tilting  term  in  the  former  resulting  in  the  production  of  turbulent  scales,  and  the 
dominance  of  turbulence  destruction  through  the  action  of  viscosity  in  the  latter.  This  conceptual 
picture  was  supported  by  the  dominance  of  the  viscous  terms  in  the  irreversible  entropy 
production. 

Visualization  of  the  layer  at  the  transitional  state  (braid  and  between  braid  planes) 
revealed  the  existence  of  convoluted  regions  of  large  density  gradient  magnitude;  no  such 
convolutions  had  been  observed  in  the  layer  that  did  not  reach  transition.  Moreover,  it  was  also 
observed  that  parcels  of  heptane  seemed  to  have  ‘broken  off  from  the  lower  stream  and  entered 
the  upper  stream.  The  striking  visual  correlation  between  the  regions  of  large  density  gradient 
magnitude,  the  desintegrated  heptane  and  the  regions  of  minimal  mass  diffusion  factor  (which 
measures  the  departure  from  mixture  ideality;  it  is  unity  for  ideal  mixtures),  provided  the 
motivation  to  examine  the  composition  and  thermodynamic  characteristics  of  the  fluid  in  the 
region  of  large  density  gradient.  Conditional  averages  revealed  that  the  fluid  in  these  regions  was 


mostly  heptane,  but  with  a  small  amount  of  nitrogen  dissolved  into  it.  Mixing  of  this  fluid  was 
observed  to  be  inhibited,  as  evidenced  by  the  smaller  that  unity  mass  diffusion  factor.  With 
increased  cutoff  for  the  proportion  of  the  value  of  the  maximum  density  gradient  magnitude  on 
w'hich  the  heptane  mass  fraction  and  mass  diffusion  factor  averages  were  conditioned,  the  fluid 
composition  became  closer  to  heptane  and  the  molecular  mixing  became  increasingly  hindered. 

Departures  from  the  perfect  gas  behavior  were  also  identified  in  contour  plots  of  the 
compression  factor  (which  is  unity  for  a  perfect  gas),  \\friereas  both  pure  streams  were  close  to 
perfect  gases,  the  mixed  fluid  showed  strong  departures  from  perfect  gas  behavior.  These 
departures  were  further  quantified  through  plots  of  the  local  difference  between  the  fluid 
temperature  and  the  local  critical  temperature,  and  the  fluid  pressure  and  local  critical  pressure  of 
the  mixture.  The  temperature  was  found  to  be  supercritical  everywhere;  therefore  the  entire  field 
is  at  supercritical  conditions  (identified  as  the  absence  of  a  two-phase  region).  However,  the 
pressure  was  found  to  be  supercritcal  in  the  pure  species  regions,  and  subcritical  in  the  fluid 
mixture  region;  this  subcritical  pressure  was  attributed  to  the  change  in  composition  of  the  fluid 
(the  critical  point  being  mixture  dependent). 

Further  studies  will  elucidate  some  aspects  of  the  assumed  PDF  representation  of  a 
supercritical  mixing  layer,  before  addressing  the  issue  of  SGS  modeling. 
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